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Abstract 

The Hubbard model on a two-leg ladder structure has been studied by a 
combination of series expansions at T = and the density-matrix renormal- 
ization group. We report results for the ground state energy Eq and spin-gap 
As at half-filling, as well as dispersion curves for one and two-hole excita- 
tions. For small U both £^0 ^-nd show a dramatic drop near t/tx_ ~ 0.5, 
which becomes more gradual for larger U . This represents a crossover from a 
"band insulator" phase to a strongly correlated spin liquid. The lowest-lying 
two-hole state rapidly becomes strongly bound as t/t_\_ increases, indicating 
the possibility that phase separation may occur. The various features are 

collected in a "phase diagram" for the model. 
PACS Indices: 71.27.+a, 71.10Fd 
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I. INTRODUCTION 



The last decade has seen a great deal of interest in spin and/or correlated electron systems 
on a ladder structure formed from two coupled chains. This work has been motivated both 
by the discovery of real materials with S = j ions forming a ladder structure and because 
ladder systems exhibit a number of interesting and surprising properties (see e.g. the review 
of Dagotto and Rice 0). 

Particularly interesting behaviour may be expected if the system includes charge degrees 
of freedom. This can be achieved by doping, to create a system of strongly correlated 
mobile holes, as in the cuprate superconductors. Systems studied include LaCu02.5 doped 
with Sr lU and Cai4Cu2404i„5 doped with Sr 0, the latter showing superconductivity under 
pressure. The t-J and Hubbard models provide alternate representations of the physics of 
such systems, and both have been studied extensively. References to most of the existing 
work on the t-J ladder system, as well as the most recent evidence for the form of the phase 
diagram can be found in Miiller and Rice 

We are here interested in the repulsive {U > 0) Hubbard model on a 2-leg ladder. 
The first study of this system, to our knowledge, was by Fabrizio, Parola and Tosatti 
who used a weak-coupling renormalization group approach to investigate the role of the 
interchain hopping in driving the system out of a Luttinger liquid phase. Earlier work, in 
this context, had considered a 2-dimensional system of weakly coupled chains P]. Further 
work 0,11 , using bosonization techniques, valid for weak interchain coupling, has identified 
a number of possible phases which the 2-leg Hubbard ladder may exhibit. In the notation 
of Balents and Fisher [§, these are denoted as CnSm where n, m represent the number of 
gapless charge and spin modes respectively. (Here n,m = 0,1,2 giving 9 possible phases). 

Numerical studies of the Hubbard ladder have been carried out for both static and dy- 
namic properties. The Density Matrix Renormalization Group (DMRG) technique has 
been used to calculate both spin and pairing correlations ||TO|,pT|, and to calculate the spin 



and charge gaps |Tl|,p!2[] as functions of the parameters of the model. The DMRG technique. 
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at least in its present form, is unable to compute the dynamic response, particulary momen- 
tum dependent properties. The one-hole spectral function has been obtained using exact 
diagonalizations (limited to 2 x 8 sites) and quantum Monte Carlo techniques . This 
latter paper also reports results for two-hole spin and charge excitations. The spin-gap at 



half-filling, obtained by DMRG on a 2 x 32 lattice |TT| showed a very sudden decrease near 
t/t±_ ~ 0.5, particularly for small U . However finite size effects are present and may mask 
the true behaviour. This work was motivated in part by a need to confirm this, and in fact 
to look for possible evidence of a phase transition near this point, as well as by a desire 
to explore the form of dispersion curves for spin and charge excitations, which cannot be 
obtained by DMRG methods. 

The Hamiltonian of the Hubbard ladder is written as 

H = -tJ2 icla,aCi+l,a,a + H.C.) - t± ^(4i,^Ci,2,a + H.C.) + U ni,a,]ni^a,i (1) 
i,a,a i,a ia 

where i labels the rungs of the ladder, a (=1,2) is a leg index, cr (=1, |) is the spin index, and 
the operators have the usual meaning. There are several instructive limiting cases. When 
U = the electrons are non-interacting and the system has two simple cosine bands of width 
At and separation 2t±. The system is metallic for all electron densities n, except for the case 
n = 1 (half-filling) and tj_ > 2t when there is a band gap and the lower (bonding) band 
is completely filled. The case is aptly referred to as a "band insulator" |0. On the other 
hand when the interchain hopping = (and U > 0) the ladder decouples to two Hubbard 
chains for which there are exact, but highly nontrivial solutions. The system is then in a 
Luttinger liquid phase, and also an insulator at half-filling. 

When both tj_ and U are non-zero there are no exact results known. However the form of 
the phase diagram can be reasonably inferred from the analytic and numerical calculations 
referred to above. A nice discussion is given by Noack et al. In the limit of large U 

doubly occupied sites are suppressed and the model reduces to a t- J ladder, with parameters 
J = 2t^/f/, J_L = 2t]_/U. The Hubbard ladder has electron hole symmetry under the 
t r ansf ormation 



Q,a,(7 1 Ci^a,a (2) 

The phase diagram and properties of the model are thus symmetric about the n = 1 case. 
We see manifestations of this in the series for various excitation energies, discussed below. 

In this paper we study the Hubbard ladder using both the method of series expansions 
]r5|^T7[ and DMRG 0. The series expansion method is complementary to other numerical 
methods and is able to provide ground state properties, excitation spectra and T = 
critical points to high accuracy. Another advantage is that one deals with a system in the 
thermodynamic limit and finite size corrections are not needed. We have used this approach 
recently in studies of the t-J model on the square lattice 0] and on a 2-leg ladder [ll9| . 



and we refer to those papers for details of the method and references to previous work. The 
emphasis of our work is on the half- filled case n = 1, and on one and two-hole excitations 
from half-filling. The series method is not well suited to handle variable electron density. 
We are not aware of any previous series work on the Hubbard ladder. Shi and Singh ||2^ 



have used a somewhat different series approach to study the Hubbard model on the square 
lattice. Their method, which introduces an artificial antiferromagnetic Ising term into the 
Hamiltonian is not appropriate here, as we do not expect any magnetic long range order in 
the ladder system. 

This paper is organized as follows. In Sec. II, we discuss briefly the methods used. In 
Sec. HI, we study the system at half-filling. In Sec. IV, we consider the system with one 
and two holes. The last section is devoted to discussion and conclusions. 

II. METHODS 

The series expansion method is based on a linked cluster formulation of standard 
Rayleigh-Schrodinger perturbation theory. We use a "rung basis" and write the Hamil- 
tonian in the form 

H = Ho + xV (3) 
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where 

Hq = -t± ^(c| 1 ,^Ci,2,a + H.c.) + UJ2 rii,a,^ni^a,i (4) 

i,a ia 

is taken to be the unperturbed Hamiltonian and 

XV = -t ^ (4a,aQ+l,a,a + H.C.) (5) 
i,a,a 

is the perturbation. Thus Hq describes decoupled rungs, and can be solved exactly while 
the intrachain hopping term couples the rungs and is treated perturbatively. 

The Hamiltonian for a single rung has 16 possible states. These are shown in Table I. 
For U < 3tj_ the lowest energy rung state is a spin-singlet containing two electrons (state 
in Table I). At f/ = 3t± there is a level crossing and the lowest energy rung state becomes 
a doublet S = ^ state with a single electron in a symmetric (bonding) state for U > 3t±. 
The eigenstates of Hq are then direct products constructed from the possible rung states. 
The ground state of Hq for the half-filling system (ra = 1) is the state that has each rung 
in a spin singlet. This is true even for U > 3tj_, as transferring an electron from a doubly 
occupied rung to another doubly occupied rung costs energy. 

To compute the perturbation series we fix the values of t± and U and expand in powers of 
X = t/t±. Without loss of generality we set t_L = 1 to define the energy scale. The series are 
then evaluated at the desired value of t using standard Fade approximants and integrated 
differential approximants pT| . 

We have also carried out large-scale DMRG calculations 0. Two different DMRG al- 
gorithms have been employed. Both are 'infinite lattice' algorithms with open boundary 
conditions. The first method uses a superblock consisting of the usual system and envi- 
ronment blocks with two added rungs in the middle. The system/environment blocks are 
augmented by one rung at a time and the superblocks always have an even number of rungs. 
The second method is similar except that only one rung is kept in the middle, meaning that 
the superblocks have an odd number of rungs. The second method allows more states to 
be retained in the blocks, whilst the even lattices dealt with in the first method are usually 
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considered to be more desirable for finite-size scaling (FSS) studies. For the first method 
typical calculations involve ladders with up to 60 rungs, keeping up to 550 states per block. 
For the second method the superb locks studied typically reached 81 rungs and up to 1500 
states were retained per block. It should be noted that the 'infinite lattice' algorithm, de- 
spite its name, can be used to obtain accurate results for finite lattices, as we shall show in 
what follows. 

III. THE HALF-FILLED CASE 

A. Ground state energy 

We first consider the ground state energy Eq at half-filling. The series method yields an 
expansion of 

Eo = Y.('n{U/t^){t/t^r (6) 

n 

The series have been computed to order for various U/tj_. The coefficients for U/t± = 8 
are given in Table II, and the other coefficients can be provided on request. The cluster data 
for this one- dimensional problem are trivial. The limiting factor is the size of the matrices 
used in obtaining the cluster energies. DMRG gives the ground state energy directly for a 
given lattice, and finite-size scaling must be used to extract the bulk result. 

In Figure ^ we show our results for the ground state energy in the half-filled case versus 
t for various U/tj_ obtained both by series expansions and DMRG calculations. For given 
t the energy increases with U as expected. For fixed U the energy decreases slowly with 
increasing t. The series are well converged for t/t_i_ ^ 0.6, but the convergence deteriorates 
rapidly at that value. The DMRG is, however, well converged even for larger t and agrees 
well with the series results for smaller t. A remarkable feature of the ground state energy is 
the sharp downturn which occurs near t/t± ~ 0.5 — 0.6 for f//t_L ^ 1, but which becomes 
smoother for larger U. For the case of free electrons {U = 0) the ground state energy is 
exactly 
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Eq 



-1, t/t^ < 1/2; 

(7) 

-1 - f [v^(2t/t^)2 - 1 - cos-i (1^)] , t/t^ > 1/2. 



and, as mentioned above, the crossover represents a transition from band insulator to a 
conductor with gapless charge and spin excitations. The quahtative similarity between the 
U = case and small U (e.g. U/t± = 0.5 in the figure) suggests that spin excitations may 
remain gapless, or nearly so, for small finite U and large t. 



We can compare our results for the ground state energy with those of Kim et al. [22 
obtained via a variational method. These are shown in Figure 2. The results are in good 
agreement for small t/t_L, but it appears that the variational method significantly over 
estimates the ground state energy for large t/t^. 

B. Spin excitations 

We now turn to the excitations, which we have computed directly using Gelfand's method 



17 1 . The lowest energy spin excitation branch results from exciting a rung from a spin 
singlet (state in Table I) to a triplet (any of states 5,6,7 in Table I). Nonzero t allows this 
to propagate coherently along the chain, giving rise to a triplet magnon excitation. Figure 
shows dispersion curves /S.s{k) for this excitation for fixed t/t±_ = 0.5 and various U/t±. 
For larger t/t± the series become too erratic. For all cases the energy minimum occurs at 
k = TT. As U increases the bandwidth of the spin excitation decreases. 

The excitation energy at A; = vr defines the spin gap = As(7r). The raw data are 
shown in Figure |^, where we plot A^/tj, versus t/t± for various U/tj_. For small t our series 



give rather precise estimates. In the limit t = the spin gap is ^iyU"^ + 16tj_ — U) which 
is 2tj_ for U = and decreases for increasing U. Near t/tj, ~ 0.4-0.5 there is a crossover 
and beyond this point the spin gap for large U exceeds that for small U. For U/t± ^ 2 this 
crossover becomes very sharp. The series extrapolations would indicate that the spin gap 
actually vanishes at t/tj_ ^ 0.5 — 0.6, but we believe this is an artefact of poor convergence 



in this region. It is expected |8;11] that at half-filling the system will have a finite spin 
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gap throughout the phase diagram for any U > 0, due to 'umklapp' processes. It is clear 
however that the spin gap becomes quite small for t/t± 0.5. 

We have also used our DMRG algorithms to calculate As. In order to obtain insight 
into the potential limitations of DMRG methods for this problem, we consider the finite-size 
scaling of in the exactly solvable U = case. As mentioned, for U = 0, the system is 
gapped for t± > 2t and gapless for t± < 2t. In Fig. § we plot A^ as a function of l/Nj-^^g for 
three values of the ratio t/t±. For t/t± = 0.49 (just below the critical ratio), the system has a 
small spin gap, and As scales smoothly to its bulk value, the finite-size corrections vanishing 
as (as is to be expected for a gapped system with open boundary conditions). For 

t/t± = 0.51 (just above the critical ratio), the gap suddenly begins to display a "sawtooth" 
dependence on A^rung, with the points of the saw scaling towards zero as A'rung oo. This 
behaviour becomes more pronounced as we move further beyond the critical ratio. It appears 
that multiple crossovers are occurring, with new states crossing over to the bottom of the 
spectrum as Arung increases. These sawtooth oscillations are presumably due to the 2-band 
structure: every so often it is is energetically favorable for an extra electron to be added 
to the top band rather than the bottom band, and a crossover occurs. This oscillatory 
behaviour hs not been exhibited before, as far as we are aware. 

In order to calculate the bulk spin gap successfully with DMRG, we must first accurately 
calculate As for a number of lattice sizes A^rung? and then extrapolate to the bulk limit, 
assuming some scaling ansatz for the finite Arung corrections. An example of this is shown 



in Table |T|, where we compare DMRG and exact results for As for finite lattices and the 
bulk limit. In this case, where the FSS is smooth, it is possible to obtain a reasonable 
estimate for the bulk spin gap, the real error being around 0.1%. Two potential problems 
emerge from our studies of gapless systems in the U = case for t/tj_ > 0.5. Firstly, because 
of the oscillating FSS it might be difficult to obtain accurate DMRG results for very large 
lattices in a regime where the gap is small or vanishing. Convergence of finite lattice results 
with the number of states, m, retained per block, must be monitored over a large range of m 
values. For a given lattice size, improved DMRG estimates can be obtained by using a finite 
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lattice algorithm However, even if highly accurate results are available for a number 
of lattice sizes, a second and more pressing problem is that, as a result of the oscillations 
or erratic FSS, extremely large lattices may be needed in order to reach a regime where a 
suitable scaling ansatz can be reliably used to extrapolate to the bulk limit, as can be seen 
for the gapless cases in Fig. ^. 

Fortunately, the presence of electron repulsion smooths out the FSS, allowing reliable 
DMRG estimates of some way beyond t/t±_ = 0.5. In Fig. ^ we show the FSS of DMRG 
estimates of in the U/t± = 1 case for various t/t±. We find that up to around t/t± = 0.55, 
the corrections to the bulk results for As scale as l/N^^^^, as in the gapped case for U = 0. 
As t/t± is increased, however, we observe initially oscillatory or erratic FSS, followed by a 
crossover to linear dependence of A^ on as might be expected for a system with 

a small gap. This is illustrated in Fig. ^ for the t/t± = 0.6 case, where a crossover to 
linear behaviour is observed as the lattice size reaches around 30 rungs. As mentioned, it 
is important to first assess the DMRG convergence for the finite lattices before attempting 
extrapolations. In Table ^V\, we show the DMRG convergence of A^ with m for the odd- 
rung algorithm, with m ranging from 200 to 1500. It can be seen in this case that the finite 
lattice estimates are sufficiently well converged to afford reliable extrapolations. In Fig. |^ 
results from the odd-rung algorithm with m = 750 are plotted along with results from the 
even-rung algorithm with m = 800. Good agreement can be seen between the two sets of 
results in the linear regime. In order to make bulk estimates for t/tj_ > 0.6, we assume a 
linear scaling ansatz. Presumably, if the resulting estimate is non-zero, as is the case, e.g., 
for t/tj_ = 0.6, there should be a crossover from linear to quadratic scaling for larger lattices 
still, in which case the estimates are (tight) lower bounds on the spin gap. For the case of 
U/tj_ = 1, we can carry this out up until t/tj_ reaches around 0.75. Beyond this value, the 
FSS is too erratic to permit accurate finite lattice estimates of the spin gap for sufficiently 
large lattices. In the inset to Fig. 5 this erratic FSS is shown for the t/t± = 0.8 case. 

Plots of the DMRG bulk estimates of A^ as functions oft/tj_ are included in Fig. ^. The 
spin gap obtained in this way is indistinguishable from the series values for small t. However, 
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the DMRG also provides well converged results for larger t provided that U is not too small. 
In the U/t± = 2 case, for example, undergoes a very rapid decrease up till t/t± ~ 0.6, 
and then flattens out at a small but finite value. The same behaviour was already seen in 



the previous DMRG calculation of Noack et al. |TT]]. Noack et al. |[TT|] did their calculation 
at a fixed lattice size 2 x 32 sites, however, whereas we performed a careful finite-size scaling 
analysis to confirm that the spin gap remains finite, rather than scaling to zero in the bulk 
limit. 

In the case U/t±_ = 1 similar effect occurs, but the FSS behaviour becomes 

completely erratic and the extrapolations fail at around t/t_i_ ^ 0.75. In Fig. |^ we depict the 
region in the {U/t±)-{t/t±) plane where the DMRG either fails to give reasonable results 
due to erratic FSS, or indicates a very small or vanishing gap, by marking it with an "F". 
Also included in Fig. |^ is the line where the series appear to give a vanishing spin gap. This 
line indicates the crossover in the physics of the system from band-insulator to strongly 
correlated Mott insulator. 



IV. 1 AND 2-HOLE STATES 

As mentioned in Section | there is considerable interest in the doped Hubbard ladder 
where the electron density n < 1. The series method is not well suited to study the effect 
of finite doping. However we are able to compute the ground state and excitation energies 
when the system contains one or two holes. 



A. 1-Hole Case 

In the t = limit a single hole will change one of the singlet rung states into an 5^ = | 
or — I bonding state, with an energy increase of —t± — Ai. Finite t will allow the hole to 
propagate along the ladder, giving a quasiparticle band. Figure |^ shows the quasiparticle 
excitation energy Aih/t_L as a function of wavenumber for the case t/t± = 0.5 and various 
U/tj_. For f/ = we have the exact results 
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Aih(fc) =tx + 2tcos(fc) (8) 

and this is seen through the vanishing of all higher terms in the series. For the choice 
t/tx = 0.5 this gives Aih = at A; = tt. We also show, for comparison, the approximate 
dispersion curve for U/t± = 4 obtained by Endres et al. through their "local rung 
approximation". This appears to overestimate the energy by about 0.4tx, although the 
overall shape is very similar. 

For increasing U the energies are depressed and there is a small decrease in the quasiparti- 
cle bandwidth although the overall cosine shape remains. The minimum of the quasiparticle 
spectrum occurs at = vr throughout. The energy zero is taken to be the ground state 
energy of the half-filled ladder. Each quasiparticle band crosses the zero level, indicating 
that the overall energy of the system is reduced when a hole is created, and this is also 
marked in Fig. |^. Figure ^ shows the "quasiparticle g function of t/t^ for various 

U/t±_ obtained from the series and from DMRG. Both methods agree up to t/t_L ~ 0.5, 
beyond which the series fail to converge. This is the same crossover region seen in the spin 
gap studies. The DMRG results suggest an upturn or change in slope for larger t/t±_ and not 
very large U , as shown in the Figure for f/ =1. However this remains somewhat speculative 
as the scaling behaviour of the DMRG estimates is unusual. We illustrate this in Figure 
^, where we plot the gap (for f//t_L = 1) versus l/N"^. The behaviour shows a qualitative 
change between t/t± = 0.5 and 0.6 becoming oscillatory and with a change of sign of limiting 



slope. Nevertheless the larger N points scale quite well. Noack et al. []Tl| saw s similar break 
in the behaviour of the charge gap at the crossover for U/t± =4. 

B. 2-Hole Case 

We would also like to explore the system doped with two holes, to see if binding occurs 
between the holes. From Table I, one can see that at zeroth order [t = 0), the energy 



gap for two holes sitting on the same rung is hiyU^ + 16t]_ — f/), which is larger than the 



gap {yW^ + 16t5_ — U) — 2t± for two holes on different rungs. Thus it is not energetically 
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favourable to have two holes on the same rung. Our present series methods, unfortunately, 
are unable to treat the latter case of two holes on separate rungs: so here we restrict ourselves 
to exploring whether the former state becomes bound at larger t. One starts from a state 
with both holes on a single rung and all other rungs in spin singlet states; and then the 
hopping term allows these holes to move along the ladder, and so generates a dispersion 
relation for this state. For example at U/t± = 8, our second order series result for the 
dispersion is 

A2h{k)/ti_ = 0.4721 + (1.805 + 1.U7 cos k){t/t±f (9) 

One can see that the minimum energy is at = vr, rather than 0, due to the fact that this is 
not the lowest energy two-hole state. Comparing with the one-hole dispersion relation, we 
find no evidence for binding of this state at higher t. 

Our series results do not preclude the existence of a two-hole bound state of more complex 
structure, and there are indications from other work that this can occur. Noack et al. |TO|JTl[| 
find, for a 2 x 32 lattice, a small binding energy {E^ ~ 0.14) for two holes, and a pair 
wavefunction which has roughly equal amplitudes on a rung and between nearest neighours 
along the same leg of the ladder. Finite size effects may be large enough to mask the true 



behaviour. Kim et al. p2[ also discuss pair formation using both DMRG and a variational 
method. They conclude that, at least for large [/, holes favour adjacent rungs to minimise 
the Coulomb energy. Our series method is unable to explore such complex pair states. 

Instead, we have used the DMRG method to compute the minimum energy of the system 
with one and two holes, up to ladders of size 2 x L {L = 60), and hence the binding energy 
defined by 

Ef, = 2[Eo{L, L - 1) - Eo{L, L)] - Eo{L - 1, L - 1) + Eo(L, L) (10) 

where Eq{N^,Ni) is the ground state energy with {Ni) up (down) electrons. We show 
this as a function of t/tx for U /t± = 2, 8 in Figure 0. As can be seen, there is no binding 
for small t/t±, and the binding energy is zero, corresponding to an unbound, well seperated 
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pair of holes. For t/t±_ larger than a critical value {t/t±)c, however, hole binding becomes 
energetically favourable. In fact as we can see from Figure 0, Ef, increases from zero 
extremely quickly once past the critical {t/t±)c. Now the binding of two holes is expected to 
be a necessary though not sufficient precondition for the phenomenon of phase separation 
(binding of many holes). Figure |TU| leads us suppose that phase separation may well occur 
for large t/t±; or else at fixed t/t± = 1, phase separation may well occur at small U/t_i_. This 
would accord with our current knowledge of the t — J ladder, where the exchange coupling 
J ~ 4t^/f/, and phase separation occurs at small U/t, i.e. small t/J. The boundary for 



two- hole binding is shown in our "phase diagram" , Figure No explicit search for phase 
separation has yet made in the Hubbard ladder, as far as we are aware. 

V. DISCUSSION AND CONCLUSIONS 

We have made the first application of the series expansion method to the Hubbard 
model on a two-leg ladder, and also used the DMRG method to explore its properties at 
half-filling at T = 0. Our series approach starts from a basis of rung states, appropriate for 
small values of the chain hopping parameter t, and we obtain perturbation series for various 
quantities in powers of t up to typically 10 terms. The series are well behaved out to typically 
t/tj_ ~ 0.6. but the convergence becomes problematical beyond that point. Our results are 
complementary to both analytical weak coupling calculations and other numerical (DMRG 
and QMC) results. 

We have calculate the ground state energy and the triplet spin excitation energies at half- 
filling, and the excitation energies of one-hole and two-hole states relative to the half-filled 



case. At half-filling the system is believed to be in a spin-gapped insulating phase for 
all values of the parameters U,t±. Our results confirm those of Noack et al. [|1^], showing 
a sharp crossover at smaller U/t_i_ from strongly correlated spin-liquid or Mott insulator (as 
for a simple Hubbard chain) to "band insulator" behaviour as t± is increased. The spin 
gap becomes very small in the spin-liquid region. We have peformed a careful finite-size 
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scaling analysis using DMRG to show that the spin gap remain finite, however, rather than 
scahng to zero in the bulk limit, at least for those moderate values of U /t±_ where a definite 
statement is possible. 

A single hole doped into the Hubbard ladder will propagate as a well defined quasipar- 
ticle, and we have compute the dispersion relation for this quasiparticle. It was found that 
the lowest energy two-hole state is not amenable to our series approach; but the DMRG 
calculations show that two holes become strongly bound at larger t/tj,. This raises the 
question of whether phase separation will occur in this region, as it does in the t — J ladder 
||23| . This question must await future investigations. 
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TABLES 



TABLE L The sixteen rung states and their energies, where m = i^l + C//^[/2 + IGt^ , 
= + Xi = l{U- ^JU^ + ml), X2 = l{U + ^JU^ + ml), and represent 



hole, I (I) represent up (down)-spin electron, J represent an electron pair. 



State label 


Eigenstate 


Eigenvalue 


Name 





^i(in)-UT))-^2(i:o)+|o|)) 


Ai 


singlet 


1 


;^(|oi)+Uo)) 


-t± 


electron- hole bonding (S'tot = ^5) 


2 


75(|0T)+|T0)) 


-t± 


electron-hole bonding {S^^^ = 5) 


3 


^(lU)- lU)) 


U-t± 


three-electron antibonding {S^ot — ~h) 


4 


TidiT)- it:)) 


u-t± 


three-electron antibonding (^tot — \) 


5 


Ui) 





triplet {Sl^ = -1) 


6 


7i(in)+ lit)) 





triplet = 0) 


7 


ITT) 





triplet {SI, = 1) 


8 


1 00) 





hole-pair singlet 


9 


73(110)- |0T)) 


u 


an electron pair and a hole singlet 


10 


III) 


2U 


an electron-pairs 


11 


^jdoD-iJ. 0)) 


l± 


(^k-ctrou-holo aniibcmdiug' (.S^^, = 


12 


;}5(|0T)-|T0)) 


t± 


electron-hole antibonding {S^^, = \) 


13 


^(IIT)+ liD) 


u + t± 


three-electron bonding {S,^, = ~i) 


14 


^(IIT)+ ITT)) 


u + t± 


three-electron bonding {S^^, = 5) 


15 


n2(|Ti)-UT)) + m(IT0)+|0l)) 


M 


mixed singlet 
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TABLE II. Nonzero coefHcients of {t/t±)^ for the ground state energy per site Eo/N, the spin 
gap As, the 1-hole energy Aih(7r), and the 2-hole energy A2h(7r) at U/t± = 8. 



n 


Eo/NJ± 




A2h/J± 


n 


Aih/J± 





-2.360679774998x10" 


1 


4.721359549996x10" 


1 


4.721359549996x10-1 





-5.278640450004x10" 


1 


2 


-8.944271909999x10- 


2 


-7.478019326001x10" 


1 


3.577708764000x10-1 


1 


-1.447213595500 




4 


-7.055810863547x10" 


2 


2.984040157758x10" 


1 


2.117784192484x10-1 


2 


-5.681892697507x10" 


2 


6 


-5.448947627757x10" 


2 


2.126553951831x10" 


1 


-1.022223221785 


3 


4.190688837075x10" 


1 


8 


-3.092587547718x10- 


2 


1.030100155673x10- 


1 


-1.054110318250x101 


4 


-2.186544916675x10- 


2 


10 


5.006769301333x10- 


3 


-2.047688788511x10- 


2 


-8.623391111290x101 


5 


7.778879476658x10" 


1 


12 


4.598641784004x10" 


2 


-2.168399847938x10" 


1 


-6. 784316243528 xlO^ 


6 


-1.450073605871 




14 


7.235038860679x10" 


2 








7 


4.979933826800 
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8 


-1.3935840289131 
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9 


4.1152389626051 





TABLE III. Comparison between exact and DMRG results for the spin gap As/t± for various 
lattice sizes Nmng for the case of U = 0, t/t± = 0.49. The DMRG results are obtained using the 
odd-rung algorithm, retaining around m = 700 states per block, and including the ground and first 
excited states in the density matrix with equal weights. The A^rung = oo DMRG result is obtained 
by extrapolating the finite Nj-^ng DMRG results over the range 63 < Nj-^^g < 81, assuming that 



the corrections to the bulk result scale as l/N^^^ 



N 

-"rung 


DMRG 


Exact 


3 


0.614071 


0.614071 


9 


0.136013 


0.135929 


19 


0.064140 


0.064131 


37 


0.046701 


0.046694 


55 


0.043092 


0.043083 


67 


0.042103 


0.042091 


81 


0.041456 


0.041438 


oo 


0.040055 


0.040000 



18 



TABLE IV. DMRG convergence of the spin gap As/t± with m, the number of states retained 
per block, for the odd-rung algorithm, where the ground and first excited states are included in 

the density matrix with equal weights. The parameters used are U/t± = 1, t/t± = 0.6. 





}11 = 


ill = /) i u 


/// = OoO 


rn = I '){) 


/// = izOU 


rn = loUli 


3 


0.2371380 


0.2371380 


0.2371380 


0.2371380 


0.2371380 


0.2371380 


7 


0.1115510 


0.1115030 


0.1115060 


0.1115040 


0.1115010 


0.1115010 


13 


0.0675853 


0.0673585 


0.0672920 


0.0672718 


0.0672626 


0.0672597 


19 


0.0475570 


0.0475439 


0.0475508 


0.0475459 


0.0475420 


0.0475383 


25 


0.0382940 


0.0383589 


0.0383977 


0.0384008 


0.0384094 


0.0384077 


31 


0.0333642 


0.0335102 


0.0336259 


0.0336437 


0.0336687 




35 


0.0314530 


0.0317155 


0.0319372 


0.0319716 


0.0320189 
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FIGURES 

FIG. 1. The ground state energy per site as function of t/t±_ for U /t±_ = 0,0.5,1,2,4,8. 
The sohd hues are extrapolations of the series using different integrated differential approximants, 
while the points connected by dashed lines are the results of DMRG calculations. Also shown are 
the exact results for {7 = 0, and the results of a variational approach [22| for U/t±_ =4,8 and 
= 0.5, 1 (full circle points). 

FIG. 2. Spin excitation spectra for t/t±_ = 0.5 and various U/t_\_, obtained from series expan- 
sions. 

FIG. 3. The spin gap versus t/t±_ for various U /t^. The dashed lines are extrapolations of 
the series using different integrated differential approximants, while the points connected by solid 
lines are the results of DMRG calculations. The filled circles (crosses) denote results obtained from 
the even-rung (odd-rung) algorithm. 

FIG. 4. Finite-size scaling of the spin gap in the noninteracting case {U = 0) for t/t±_ = 0.49 
(a), t/t_\_ = 0.51 (b) and t/t_i_ = 0.55 (c). These are exact results for open lattices. 

FIG. 5. Finite-size scaling of the spin gap A^ for U/t±_ = 1 and t/t_\_ = 0.54,0.6, obtained 
from DMRG calculations. For t/t±_ = 0.6, the results from the two different DMRG algorithms: 
odd-rungs (crosses) and even-rungs (full dots), retaining 750 and 800 states per block respectively, 
are compared. The inset shows the FSS in the U/t±_ = 1, t/t±_ = 0.8 case using the odd-rung 
DMRG algorithm with m = 1200. 

FIG. 6. Phases and critical lines for the Hubbard ladder at half-filling, in the plane of U /t_\_ 
versus t/t±_ (see text). The phase boundary between spin liquid and band insulator is obtained 
approximately by determining the position where the spin gap vanishes according to the series. 
The region that two holes bind (do not bind) is marked by B(NB). The region in which DMRG 
fails to determine whether the system has a spin gap, due the irregular finite size scaling, is marked 
by F. Also marked is the region that the 1-hole gap is positive/negative. 
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FIG. 7. Excitation spectra for one-hole bonding states for t/t±_ = 0.5 and various U/t_\_, 
obtained from series expansions. Also shown are the results of the local rung approximation ||l^ 
at U /t±_ =4 (dashed line) . 

FIG. 8. 1-hole "quasiparticle gap" as a function of t/t±_ for various U /t±. The dashed lines are 
the extrapolations of the series using the different integrated differential approximants, while the 
points connected by solid lines are the results of DMRG calculations. Also shown are the exact 
results for the case C/ = 0. 

FIG. 9. 1-hole "quasiparticle gap" Ai/,(7r) for U/t^ = 1 and t/t^ = 0.5,0.52,0.54,0.6, as a 
function of l/A^^^ungi obtained from DMRG calculations. 

FIG. 10. Binding energy of two holes versus t/t±_ for various U/t±, obtained from DMRG 
calculations. 
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